##########################################################################################

library(dplyr)
library(data.table)
library(optparse)
library(ggplot2)
library(ggsci)
library(ggrepel)
library(ggpubr)
library(parallel)

##########################################################################################
option_list <- list(
    make_option(c("--sample_list_file"), type = "character"),
    make_option(c("--input_file"), type = "character"),
    make_option(c("--out_path"), type = "character")
)

if(1!=1){
    
    ## 突变信号文件
    sample_list_file <- "~/20220915_gastric_multiple/rna_combine/analysis/config/tumor_normal.list"

    ## 免疫浸润
    input_file <- "~/20220915_gastric_multiple/rna_combine/analysis/images/ddr_gsva/gsva_ddr.tsv"    

    ## 输出
    out_path <- "~/20220915_gastric_multiple/rna_combine/analysis/images/ddr_gsva"

}

###########################################################################################
parseobj <- OptionParser(option_list=option_list, usage = "usage: Rscript %prog [options]")
opt <- parse_args(parseobj)
print(opt)

sample_list_file <- opt$sample_list_file
input_file <- opt$input_file
out_path <- opt$out_path

##########################################################################################

info <- data.frame(fread(sample_list_file))
dat_cibersort <- data.frame(fread(input_file))

##########################################################################################

class_type <- c( "Normal" , "IM" , "IGC" , "DGC")

rownames(dat_cibersort) <- dat_cibersort$Sample
dat_cibersort <- t(dat_cibersort[,2:ncol(dat_cibersort)])

##########################################################################################
## 按人的不同病理类型合并样本
## 若一个人同一病理类型多个样本，均中位数
dat_tpm_all <- Reduce(function(x,y)bind_cols( x , y),mclapply(unique(info$ID) , function(id){

    tmp_info <- subset( info , ID == id )

    result <- data.frame()
    ## 若一个人同一病理类型多个样本，均中位数
    for(class in unique(tmp_info$Class)){
        tmp_sample <- tmp_info[tmp_info$Class==class,"Tumor"]
        tmp_tpm <- dat_cibersort[,tmp_sample]

        if(length(tmp_tpm) > 22){
            value <- apply( tmp_tpm , 1 , median )
        }else{
            value <- tmp_tpm
        }
        
        result_tmp <- data.frame( sample = value )
        colnames(result_tmp)[1] <- paste0( id , "_" , class)

        if( nrow(result) > 0){
            result <- cbind(result , result_tmp)
        }else{
            result <- result_tmp
        }
    }

    ## Normal样本
    tmp_sample <- unique(tmp_info$Normal)
    if(tmp_sample!="#N/A"){
        tmp_tpm <- dat_cibersort[,tmp_sample]
        value <- tmp_tpm
        result_tmp <- data.frame( sample = value )
        colnames(result_tmp)[1] <- paste0( id , "_Normal")

        ## 输出结果
        result <- cbind(result , result_tmp)
    }
    
    result

},mc.cores=10))

## 只关注Normal、IM、IGC、DGC的表达情况
col_names <- grep( paste( class_type , collapse="|") , colnames(dat_tpm_all) , value = T )
dat_tpm_all <- dat_tpm_all[,col_names]

##########################################################################################

igc_class <- c("IM + IGC + DGC" , "IM + IGC")
dgc_class <- c("IM + IGC + DGC" , "IM + DGC")

im_igc_sample <- unique(info[info$Type %in% igc_class , "ID"])
im_dgc_sample <- unique(info[info$Type %in% dgc_class , "ID"])

##########################################################################################
## 分Normal->IM->IGC，Normal->IM->DGC

sample <- sapply( strsplit(colnames(dat_tpm_all)[1:ncol(dat_tpm_all)] , "_") , "[" , 1)
class <- sapply( strsplit(colnames(dat_tpm_all)[1:ncol(dat_tpm_all)] , "_") , "[" , 2)

result <- c()
for( cell_type in rownames(dat_tpm_all) ){
    print(cell_type)

    tmp_ratio <- dat_tpm_all[cell_type,]
    ratio <- as.numeric(tmp_ratio)

    tmp_dat <- data.frame( Sample = sample , Class = class , Ratio = ratio , CellType = cell_type )

    ## Normal->IM->IGC看基因的表达变化情况
    tmp_dat_igc <- subset( tmp_dat , Sample %in% im_igc_sample )
    tmp_dat_igc <- subset( tmp_dat_igc , Class %in% c( "Normal" , "IM" , "IGC" ))
    tmp_dat_igc$Class <- factor( tmp_dat_igc$Class , levels = c( "Normal" , "IM" , "IGC" ) , order = T )
    tmp_dat_igc$Type <- "IM + IGC"

    ## Normal->IM->DGC看基因的表达变化情况
    tmp_dat_dgc <- subset( tmp_dat , Sample %in% im_dgc_sample )
    tmp_dat_dgc <- subset( tmp_dat_dgc , Class %in% c( "Normal" , "IM" , "DGC" ))
    tmp_dat_dgc$Class <- factor( tmp_dat_dgc$Class , levels = c( "Normal" , "IM" , "DGC" ) , order = T )
    tmp_dat_dgc$Type <- "IM + DGC"

    ## Normal->IM->GC看基因的表达变化情况
    tmp_dat_gc <- tmp_dat
    tmp_dat_gc$Class <- ifelse( tmp_dat_gc$Class %in% c("IGC" , "DGC") , "GC" , tmp_dat_gc$Class )
    tmp_dat_gc$Class <- factor( tmp_dat_gc$Class , levels = c( "Normal" , "IM" , "GC" ) , order = T )
    tmp_dat_gc$Type <- "IM + GC"
    ## 一个既有IGC又有DGC，合并
    tmp_dat_gc <- tmp_dat_gc %>%
    group_by(Sample , Class) %>%
    summarize( Ratio = median(Ratio) , CellType = unique(CellType) , Type = unique(Type) )

    tmp_dat <- rbind( tmp_dat_igc , tmp_dat_dgc , tmp_dat_gc )

    result <- rbind( result , tmp_dat )
}

##########################################################################################
out_name <- paste0( out_path , "/gsva_ddr.MutipleStage.pdf" )

my_comparisons_1 <- list( 
    c(1, 2), c(1, 3) , 
    c(2, 3) )

dat_plot <- result
dat_plot$Type <- factor( dat_plot$Type , levels = c("IM + GC" , "IM + IGC" , "IM + DGC") )

#dat_plot <- subset( result , CellType %in% show_cellType  )
#dat_plot$CellType <- factor( dat_plot$CellType , levels = order_cellType , order = T )

plot <- ggplot( dat_plot , aes( x = Class , y = Ratio , color = Class ) ) +
    geom_line( aes( group = Sample ) , size = 0.4 , color = "gray" ) +  ## 配对样本加线
    geom_boxplot(alpha =1 , outlier.size=0 , size = 0.9 , width = 0.6) +
    geom_jitter(position = position_jitter(0.17) , size = 1 , alpha = 0.7) +
    facet_grid(CellType~Type,space='free_x',scales='free_x') +
    scale_fill_npg()+
    scale_color_npg()+
    xlab(NULL) +
    ylab("GSVA Score")+
    theme_bw() +
    stat_compare_means(comparisons = my_comparisons_1 , label.y = c(1,1.2,1.4)) +
    theme(panel.background = element_blank(),#设置背影为白色#清除网格线
        legend.position ='left',
        legend.title = element_blank() ,
        panel.grid.major=element_line(colour=NA),
        plot.title = element_text(size = 8,color="black",face='bold'),
        legend.text = element_text(size = 10,color="black",face='bold'),
        axis.text.y = element_text(size = 10,color="black",face='bold'),
        axis.title.x = element_text(size = 10,color="black",face='bold'),
        axis.title.y = element_text(size = 10,color="black",face='bold'),
        axis.ticks.x = element_blank(),
        axis.text.x =  ,
        axis.line = element_line(size = 0.5)) 
ggsave(file=out_name,plot=plot,width=7,height=20)

out_name <- paste0( out_path , "/gsva_ddr.MutipleStage.tsv" )
write.table( dat_plot , out_name , row.names = F , quote = F , sep = "\t" )


##########################################################################################
my_comparisons_1 <- list( 
    c(1, 2) , c(1, 3) , c(1, 4) , 
    c(2, 3) , c(2, 4) ,
    c(3, 4)  
    )
tmp_dat <- unique(dat_plot[,1:4])
tmp_dat <- subset(tmp_dat , Class!="GC")
tmp_dat$Class <- factor( tmp_dat$Class , levels = c( "Normal" , "IM" , "IGC" , "DGC" ) , order = T )

out_name <- paste0( out_path , "/gsva_ddr.MutipleStage.oneImage.pdf" )

plot <- ggplot( tmp_dat , aes( x = Class , y = Ratio , color = Class ) ) +
    geom_line( aes( group = Sample ) , size = 0.4 , color = "gray" ) +  ## 配对样本加线
    geom_boxplot(alpha =1 , outlier.size=0 , size = 0.9 , width = 0.6) +
    geom_jitter(position = position_jitter(0.17) , size = 1 , alpha = 0.7) +
    facet_grid(CellType~.,space='free_x',scales='free_x') +
    scale_fill_npg()+
    scale_color_npg()+
    xlab(NULL) +
    ylab("Stem score")+
    theme_bw() +
    stat_compare_means(comparisons = my_comparisons_1) +
    theme(panel.background = element_blank(),#设置背影为白色#清除网格线
        legend.position ='none',
        legend.title = element_blank() ,
        panel.grid.major=element_line(colour=NA),
        plot.title = element_text(size = 8,color="black",face='bold'),
        legend.text = element_text(size = 10,color="black",face='bold'),
        axis.text.y = element_text(size = 10,color="black",face='bold'),
        axis.title.x = element_text(size = 10,color="black",face='bold'),
        axis.title.y = element_text(size = 10,color="black",face='bold'),
        axis.ticks.x = element_blank(),
        axis.text.x = element_text(size = 10,color="black",face='bold') ,
        axis.line = element_line(size = 0.5)) 
ggsave(file=out_name,plot=plot,width=4,height=23)

out_name <- paste0( out_path , "/gsva_ddr.MutipleStage.oneImage.tsv" )
write.table( dat_plot , out_name , row.names = F , sep = "\t" , quote = F )